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Abstract. We study thermodynamic properties and the electrical conductivity 
of dense hydrogen and deuterium using three methods: classical reactive Monte 
Carlo (REMC), direct path integral Monte Carlo (PIMC) and a quantum 
dynamics method in the Wigner representation of quantum mechanics. We 
report the calculation of the deuterium compression quasi-isentrope in good 
agreement with experiments. We also solve the Wigner-Liouville equation of 
dense degenerate hydrogen calculating the initial equilibrium state by the PIMC 
method. The obtained particle trajectories determine the momentum-momentum 
correlation functions and the electrical conductivity and are compared with 
available theories and simulations. 
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1. Introduction 

During the last decades significant efforts have been made to investigate the 
thermophysical properties of dense plasmas. The importance of this activity is 
mainly connected with the creation of new experimental facilities. Powerful current 
generators and ultrashort lasers are used for production of very high pressures 
which cannot be reached in traditional explosive devices and light-gas guns. Such 
experiments give valuable information about various properties of strongly coupled 
plasmas. This allows one to obtain a wide-range equation of state and to verify various 
theoretical approaches and numerical methods. On the other hand these results are of 
fundamental interest also for various astrophysics and solid state physics applications. 

Here we report new results on thermodynamic properties and electrical 
conductivity of dense hydrogen and deuterium. Thermodynamics is calculated by 
two approaches: the REMC method [1,2] and the PIMC approach [3]. Main attention 
is paid to the region of hypothetical plasma phase transition (PPT). We also make 
a comparison with recent experimental results on the quasi-isentropic compression of 
deuterium. 

To calculate the electrical conductivity we use quantum dynamics in the Wigner 
representation of quantum mechanics. The Wigner-Liouville equation is solved 
by a combination of molecular dynamics (MD) and MC methods. The initial 
conditions are obtained using PIMC which yields thermodynamic quantities such 



Thermodynamic properties of plasma media 



2 



as the internal energy, pressure and pair distribution functions in a wide range 
of density and temperature. To study the influence of the Coulomb interaction 
on the dynamic properties of dense plasmas we apply the quantum dynamics in 
the canonical ensemble at finite temperature and compute temporal momentum- 
momentum correlation functions and their frequency-domain Fourier transforms. We 
discovered that these quantities strongly depend on the plasma coupling parameter. 
For low density and high temperature our numerical results agree well with the 
Drude approximation and Silin's formula [4], but with increasing coupling parameter 
deviations grow. 

2. Simulation methods 

A complete description of the REMC method can be found in Refs. [1,2]. Here 
we consider only molecular hydrogen dissociation and recombination: H2 2H. 
Ionization can be neglected at temperatures lower than the hydrogen ground state 
energy (oc 13.6 eV) and at moderate densities. The effective pair potentials between 
different particle species are approximated by Buckingham-EXP6 potentials, corrected 
at small distances [5] . Our REMC simulations have been performed in the canonical 
ensemble for hydrogen and deuterium. We use 3 types of MC moves: particle 
displacement, molecular dissociation into atoms and recombination to a molecule. The 
expressions for probabilities of the two last moves are given by the internal partition 
functions of atoms and molecules Z^ 1 *, Z^f. All electrons (in atoms and molecules) are 
assumed to be in the ground state. Further, Z™ 1 contains only translational degrees 
of freedom, Zffl contains, in addition rotational and vibrational degrees of freedom. 
For the latter we numerically solve the Schrodinger equation in the central-symmetric 
field, as described in Ref. [6] , which yields the energy levels E n i . 

The PIMC method allows for first-principle simulations of dense plasmas at 
arbitrary coupling and up to moderate degeneracy parameters, for details on our 
method, see [3, 15]. It has been used for calculation of thermodynamic properties of 
hydrogen and hydrogen-helium plasmas and electron-hole plasma in semiconductors. 

Finally, we briefly describe our quantum dynamics (QMD) simulations of the 
conductivity. Our starting point is the canonical ensemble-averaged time correlation 
function [7] 

C FA {t) = /#(0)i(i)) = Z~ lr Fr ^pe i6t *' h Ae- i6t ^ n \ , (1) 

where F and A are operators of arbitrary observables and H is the Hamiltonian of the 
system which is the sum of the kinetic K and the potential U energy operators and 

the complex time is t c — t — ihf3/2, and /3 = l/ksT. Z = Tr |e~^ j is the partition 

function. The Wigner representation of ([1]) in a w-dimensional space is 

C FA (t) = {2nh)- 2v x J f dtHd^F^A^WimiMtiihp), 

where \ii = {pi,qi), {i = 1,2), and p and q comprise the momenta and coordinates, 
respectively, of all particles. A{n) and F(/i) denote Weyl's symbols of the operators 

A{ P )= J dfe-** (q-\ 
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and W{jjb\; ji2', t\ ih(3) is the spectral density expressed as 

1 
Z 



e iHt*/U 



e -iHt c /h 



91 T 



(2) 



As has been proved in [8,9], VF obeys the following integral equation: 

W{ni;nv;t;ih(i) = W(p ,qo;po,q ;ifif3) + 

\ J dr J dsW(p T - s,q T ;p T ,q T T;ih(3)w(s,q T ) - 

\ J dr J dsW (pT,q~T;PT ~ s,q T ;T;in/3)m(s,q T ), 

where W(fii; fi2] ih(3) = W{nx\^i2]Q',iTip) is the initial condition equation, which can 
be presented in the form of a finite difference approximation of the Feynman path 
integral [8,9]. The expression for W has to be generalized to account for the spin 
effects. This gives rise to an additional spin part of the initial density matrix, e.g. [15]. 
Also, to improve the simulation accuracy the pair interactions U a b, are replaced by 
an effective quantum potential Uf?, such as the Kelbg potential [10]. For details we 
refer to Refs. [3, 12-15], for recent applications of the PIMC approach to correlated 
Coulomb systems, cf. [16-19]. 

The solution of the integral equation ([3]) can be represented by an iteration series 

W* = W l + K l T W T = W l + K l Tl W Tl + . . . , 

where W* and W T1 are the initial quantum spectral densities evolving classically 
during time intervals [0,t] and [0, t{\, respectively, whereas Krl +1 are operators that 
govern the propagation from time Tj to rj+i, see e.g. [11]. Thus the time correlation 
function becomes 

c FA (t) - (0|w*) = {0r*) + {4>\k 1 ti w ti ) + ... 

where </>(/ii; ^2) = F(/ii)A(fi2) and the parentheses (...[...) denote integration over 
the phase space {/ii; ^2}- 

The iteration series for Cf_a(£) can be efficiently computed using MC methods. 
We have developed a MC scheme which provides domain sampling of the terms giving 
the main contribution to the iteration series, cf. [8,9]. For simplicity, in this work, 
we take into account only the first term of iteration series, which is related to the 
propagation of the initial quantum distribution W according to the Hamiltonian 
equations of motion. This term, however, does not describe pure classical dynamics 
but accounts for quantum effects [20] and, in fact, contains arbitrarily high powers of 
the Planck's constant. 



3. Numerical results 

3.1. Deuterium compression isentrope 

To calculate an isentrope one has to determine the entropy which is defined by: 
S = -[dF(T,V,N)/dT] v , where F = E^i-^Mi - PV is the free energy of a 
two-component system of atoms and molecules. The chemical potentials \ii of each 
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component are obtained with the the test particle method [22]: fi = p, ld + p, r , where 
p ld is the ideal gas chemical potential: p ld — k B T\og(A 3 N/V) and A is the de Broglie 
wavelength of a particle. The residual chemical potential p, r can be evaluated as [23]: 

// = -k B Tlog[L a exp(-AU/k B T)}. 

Here AU denotes the change of configurational energy produced by the insertion of one 
additional particle and L a is the ratio of allowed (nonoverlapped) insertion intervals 
along trajectories which traverse (parallel to any of the axes) the simulation box from 
side to side, to the length of the box [23] . The test particle is then inserted randomly 
into some point in the allowed intervals and the change in configurational energy AU 
is evaluated. The main advantage of this approach is the possibility of calculation of 
chemical potentials at high densities, where the usual test particle method tends to 
fail. The chemical potentials are calculated separately for atoms and molecules. 

The isentrope can also be calculated by using Zel'dovich's method [24,25]. From 
the first law of thermodynamics the characteristic equation for the temperature along 
the isentrope can be derived: 

dT _ T 

~d~V ~ ~(dE/dP) v ' 

We integrate this equation with the initial condition corresponding to an experimental 
point at low pressure taking the temperature at this point to be close to that from 
the Widom's test particle method. The coefficient (dE/dP)v is obtained from the 
interpolation functions E(T, V) and P(T, V), which are constructed from the REMC 
calculation on the grid of isotherms and isochores covering the experimental isentrope. 

Calculations were performed in a cubic simulation box with periodic boundary 
conditions, and with a cutoff radius equal to half of the box length. The initial particle 
configuration was an fee lattice for every input density p = Numn/V + Nn 2 mn 2 /V 
with = Nn 2 — 250. The system was equilibrated for 2-10 7 steps, and additional 10 7 
steps were used for the calculation of thermodynamic values. Averaging of 20 blocks 
was used to calculate the statistical error, which did not exceeded 2% for pressure and 
energy. 

The results for three shock Hugoniots of gaseous deuterium with three different 
initial densities (po = 0.1335 g/cm 3 , Pq = 1.57 GPa; po — 0.153 g/cm 3 , Pq = 
2.03 GPa; p = 0.171 g/cm 3 , P = 2.72 GPa) obtained by REMC and DPIMC 
methods can be found elsewhere [26]. Here we present our recent results concerning 
the calculation of deuterium quasi-isentrope of compression. A dramatic increase of 
the conductivity of deuterium by 4-5 orders of magnitude [27] at pressures ~1 Mbar 
and densities about 1 g/cm 3 indicates that one might also expect peculiarities in 
the thermodynamic properties. Indeed, there are experimental results on quasi- 
isentropic compression of deuterium in a cylindrical explosive chamber which show 
a 30% density jump at a pressure of about 1.4 Mbar [28]. Using the deuterium free 
energy from REMC and Widom's test particle methods we calculated the compression 
isentrope of deuterium. The results are shown in Fig. [T] and compared to experimental 
data [28] . The excellent agreement proves that no special assumptions about a phase 
transition are needed to explain such a remarkable behavior of the experimental points, 
apparantly dissociation effects mostly contribute to this phenomenon. 

A comparison of our results with other theoretical and computational methods 
is shown in the phase diagram Figj2] The temperature and density jump along the 
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isentrope is close to the predicted boundary of the deuterium phase transition from 
the molecular to the atomic phase [29-31]. The rectangle shows the region of slow 
convergence of the simulations to the equilibrium state which might be an indication 
of a thermodynamic instability. Our preliminary analysis shows that this state can 
be stable (up to ~ 2 ■ 10 9 Monte-Carlo steps are required), but to investigate the 
possibility of phase separation at these conditions one needs to apply special Monte 
Carlo methods. 
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Figure 1. Deuterium quasi-isentropic compression. Experiment: 1 — [28]. This 
work: 2 — Widom's test particle method, 3 — Zel'dovich's method. Vertical 
line shows the density at which a sharp electircal conductivity rise is observed 
experimentally [27,32], 

Figure 2. Deuterium phase diagram. Melting curve: 1 — experiment [33]; 2 — 
extrapolation of experiment [33]. Theoretical phase diagram: 3 — [34]; 4 — [35]. 
Boundary of the possible phase transition: 5 — [29]. Atomic fluid: 6 — [30,31]. 
Molecular liquid: 7 — [30, 31] . This work: 8 — isentrope; the black rectangle 
shows the region of slow convergence to the equilibrium state. 



3.2. Quantum momentum-momentum correlation functions 

We now compute the dynamic conductivity of a strongly coupled hydrogen plasma. 
The results obtained were practically insensitive to the variation of the whole number 
of particles in the Monte Carlo cell from 30 up to 60 and also to the number n of 
high temperature density matrices in the path integral representation of the initial 
state which ranged from n = 20 to 40. Estimates of the average statistical error 
gave a value of the order 5-7%. According to the Kubo formula [7] our calculations 
include two stages: (i) generation of the initial conditions (configuration of protons 
and electrons) in the canonical ensemble with the probability being proportional to 
the quantum density matrix and (ii), generation of the dynamic trajectories in phase 
space, starting from these initial configurations. 

First we discuss the momentum- momentum autocorrelation functions (MMCF) 
which are shown in Fig. [3] for various temperatures and densities. The plasma density 
is characterized by the Brueckner parameter r s defined as the ratio of the mean 
interparticle distance d — ( 47r ( ra 3 +ra \ ) 1 ^ 3 to the Bohr radius as, where n e and n p 
are the electron and proton densities. The monotonic decay of the MMCF at low 
density transforms into aperiodic oscillations at high densities. The tails of the MMCF 
clearly show collective plasma oscillations. The damping time of the initial decay of 
the MMCF is strongly affected by variations of density and temperature. At constant 
density the decay time is at least two times smaller for T = 200 000 K compared to 
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T = 50 000 K. As it follows from Fig. [3] the damping time is also sensitive to the 
density. The damping time of the MMCF increases when the density increases. The 
physical reason is the tendency towards ordering of the charges. 




Figure 3. Typical MMCF in canonical ensemble for different densities (r s ) 
and three temperatures: T = 200 000 K (left), T = 100 000 K (central) and 
T = 50 000 K (right). Time is presented in atomic units. 



3.3. Electrical conductivity 

Figure |4] presents the real part of the diagonal elements of the electrical conductivity 
tensor versus frequency computed from the real part of the Fourier transform of 
the MMCF which characterizes the Ohmic absorption of electromagnetic energy. 
Collective plasma oscillations give the main contribution in the region of lo/lo p ~ 1, 
where w 2 = 4irn e e 2 /m e is the plasma frequency. Reliable numerical data in this 
region require very long dynamic trajectories. Their initial parts are presented in 
Fig. [3] The high frequency tails of the dynamic conductivities coincide with analytical 
Drudc-like expressions for fully ionized hydrogen plasma obtained in [4, 36]. For 
low frequency analytical estimations are going to infinity and this is the reason of 
discrepancy between numerical results and analytical estimations. With increasing 
plasma density non-monotonic behavior of the dynamic conductivity in the region of 
several plasma frequencies is observed. These oscillations can be an indication of the 
transparency window (low absorption coefficient) of the strongly coupled hydrogen 
plasma. 

The agreement with Drudc-like formulas for weakly coupled plasma is due to 
the fact that the main contribution to the high frequency region comes from the 
fast trajectories with high (virtual) energy. This means that interaction of electrons 
with each other and protons only weakly disturbs the behavior of the high-frequency 
tails of dynamic conductivity in comparison with ideal plasma. Now let us consider 
the dynamic conductivity at very low plasma density, namely when the Brueckner 
parameter r s is equal to 43.2. Fig. \5\ presents MMCF and dynamic conductivity in a 
wide region of temperatures from T = 50 000 K up to T — 10 000 K. At temperatures 
lower than T = 50 000 K the plasma consists mainly of atoms. As a consequence 
the initial fast decay of the MMCF is modulated by the high frequency oscillations 
related to the motion of electrons inside the atoms. So the high frequency tail of the 
dynamic conductivity has a new maximum associated with the condition Ry — Tko. 
As it follows from Fig. [5] the position of this peak strongly depends on temperature 
and is shifted to lower frequency when temperature decreases. The physical reason 
is the growth of the energy levels population of hydrogen atoms with the increase of 
temperature. Analytical estimations for fully ionized plasma which are also presented 
in Fig. [5] give essentially larger values of dynamic conductivity. 
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Figure 4. Real part of the Fourier transform of the MMCF (line 1) for densities 
related to three temperatures: T = 200 000 K (left column), T = 100 000 K 
(central column) and T = 50 000 K (right column) and to r s = 30 (top row) 
and r B = 10 (bottom row). Points 2 and 3 present analytical results obtained 
according to [4, 37] respectively. Frequency and scaled dynamic conductivity is 
given in units of plasma frequency. 




Figure 5. Typical MMCF in canonical ensemble for r s = 43.2 and four 
temperatures: T = 50 000 K, T = 20 000 K, T = 15 000 K and T = 10 000 K, 
time is presented in atomic units (left figure). Real part of the Fourier transform 
of MMCF (line 1) for density related to r s = 43.2 and temperatures T = 10 000 K 
(central figure) and T = 20 000 K (right figure). Data points 2-5 present analytical 
results obtained according to [4,36-38] respectively. Arrow relates to the condition 
Ry = hui. 



4. Conclusion 

We have used three approaches for the investigation of thermodynamic properties and 
electrical conductivity of dense hydrogen and deuterium plasma. Using two different 
methods of isentrope reconstruction from simulation results we obtained very good 
agreement with experimental data on quasi-isentropic compression of deuterium. We 
also applied the quantum dynamics approach to hydrogen plasma in a wide region 
of density and temperature. Calculating the MMCF we determined the dynamic 
electrical conductivity and compared the results with available theories. Our results 
show a strong dependence on the plasma coupling parameter. For low density and 
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high temperature the numerical results agree well with the Drude approximation, 
while at higher values of the coupling parameter we observe a strong deviation of 
the frequency dependent conductivity and permittivity from low density and high 
temperature approximations. 
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